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ABSTRACT 

A   computer   model  is   constructed  for   the  purpose  of 
simulating  interception  of   a   near    earth   satellite  by  a 
ground   launched  ballistic   interceptor   with   emphasis   on 
adaptability   of   the   model   to   the  investigation  of   interceptor 
performance  in  terms   of  interceptor  and   satellite   character- 
istics.     The  model  is  keyed  to   instantaneous   satellite 
positions,    and   provides   on   call  information  about    minimum 
interceptor   velocity  requirements,    parameters    of   possible 
interceptor  trajectories  and  relative  velocities  at  inter- 
ception for  those  trajectories.      The  model  is   constructed 
with  sufficient  generality  and  flexibility  so  that  it   may  be 
readily  adapted  to   simulation  or  computation  of   other 
intercept   problems. 
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CHAPTER   I 
INTRODUCTION 

The  fundamental    equations   of   motion  for   a   body- 
acting   under   the   influence   of   a   central   force  have  been 
well   known  and   exhaustively   treated   since   the   times   of 
Newton  and   Kepler.      These   equations,    which   will   not  be 
red&rived  here,    form   the  basis   of   the   model   to   be   devel- 
oped.     It   should  be  noted,    however,    that   this    simple 
formulation  involves    some  assumptions   that   do   not   com- 
pletely  correspond   to    the  physical   environment   of   a   near 
earth    satellite.       These  assumptions   and   their    effects    on 
the   model   must   be   considered. 

BASIC  ASSUMPTIONS 

1.  The   earth  and   the   satellite  are   considered  as   point 
masses    with   the   mass   of   the   satellite  negligible   compared 
to   the   mass    of   the   earth. 

2.  The   earth's   gravitation  is   the  only  force   acting   on 
the   satellite. 

3.  The   coordinate   system    is   fixed  in  space. 

DISCUSSION  OF  ASSUMPTIONS 

If   the   earth   were  perfectly   spherical   with   its    mass 


distributed    symmetrically  about  its    center,    all  "its   mass 
could   be   considered   to   be   concentrated  at    one   point. 
However   the   earth   is    slightly   oblate   to   a   first   approx- 
imation and   there  are   slight     asymmetries   in  the   distri- 
bution of   its    mass,    particularly  in  the   surface   distribution 
of    continental   and   ocean   masses.       This   oblateness   is 
responsible   for   the   major   anomalies   in  its   gravitational 
field,    which   affect   the  behavior   of   near   earth   satellites 
in  two   important   ways  [3J  .      First,    it   causes   the  orbital 
plane  to   rotate  about   the   earth's   axis   in  the  direction 
opposite   to   the   satellite's    motion.      This   perturbation  is 
known  as   regression  of   the  nodes   and  its   rate   is   given  by 

M±=   -9.97  £L|      "      (l-e2)-2   cos   i   degrees/day  (1.1) 

where   R   is   the  radius   of   the   earth,    a   is   the   semi-major 
diameter   of   the   orbit,    e  is   the   eccentricity  of   the  orbit 
and   i   is   the   inclination   of   the   orbital   plane   to   the   plane 
of   the   earth's    equator.       These   terms   and   symbols    are 
discussed  in   detail   in   Chapter   II.      The   notation   through- 
out  this   paper   will   be  in  accordance   with   the   table   of 
symbols   unless   otherwise   specified.       Second,    the   earth's 
oblateness    causes    the   perigee   position  of   the   orbit   to 


rotate  about   the  principal   focus   of   the   orbit   with   the 
rate  being   given  by 

(ion;  r>        ^  p 

Rj     '       (l^e    )~      (5    cos   i-1)    degrees/day.  (1.2) 
a  I 

The  atmosphere  offers   increasing  resistance   to   the 
motion  of   a   satellite  as   its   altitude   decreases.      The 
primary   effect   of   atmospheric   drag  is   to   reduce  the   major 
diameter   of   the   orbit,    mostly  at   the   expense  of   apogee 
distance.      This   makes   the   orbit    more  nearly   circular   with 
the  passage   of   time   while  reducing  the   satellite's   period 
of   revolution  and  increasing  its   velocity  in  accordance   with 
Kepler's   third  law. 

In  addition  to   these   major    effects,    satellite   orbits 
are  perturbed  to  a   much   smaller  degree  by  solar   wind, 
radiation  pressure  and  the   gravitational  attraction  of 
other   celestial  bodies,    especially   the   moon. 

For    most   applications    of   the   model   to   be  developed, 
the   errors   mentioned  above   may  be   disregarded.       The 
model   is   geared  to   inspection  of   instantaneous    satellite 
position  and   velocity   vectors,    and  orbital   parameters    may 
be  considered  to   be  no   older   than  one  period   of   revolution. 
However,    the   model   should  not   be  applied  to    orbits  having 


altitudes   of   less   than   100    miles   because   of   the  increased 
effects   of  atmospheric   drag,    nor   to   orbits    with   apogee 
distances    of   several   earth   radii  because   of   perturbations 
due  to   the   moon's   gravitational  influence.      Neither   should 
the  model  be  used  to   predict   the  position  of   the  satellite 
for   more   than  one   or   two   orbital  periods  because   of   the 
cumulative   errors   indicated  by   equations    (1.1)    and    (1.2). 
These  restrictions   are  not    severe,    for   in  fact   most 
objects    of   interest   in  the  interception  problem,    such  as 
reconnaissance  and  weapon  carrying  satellites   will  orbit  at 
altitudes   between   100   and   1000    miles   due  to   the  require- 
ments  of   their   missions.      Substitution  of  these  values 
into   equations    (1.1)  and  (1.2)    reveal  that   even  in  the 
worst   cases,    the  perturbations  amount   to  only  a  fraction 
of  a   degree  per  revolution. 


CHAPTER    II 
ORBITAL    MECHANICS 

Kepler's   brilliant  formulation  of  the  laws  of   planetary- 
motion   early  in  the   seventeenth   century   were  provided 
with   a   sound   theoretical  basis   by   Newton's   laws   of    motion 
and   were   shown  to   apply   generally   to   bodies   acting  under 
the  influence  of  a   central  force.      These  laws   can  be 
stated  for   an   earth   satellite  as   follows: 

1.  The   orbit   of   a   satellite  is   an   ellipse   with   the 
center  of  the  earth  at  one  of  its  foci. 

2.  The  radius   vector   of   a   satellite   describes    equal 
areas   in   equal   times. 

3.  The   square  of   the  period   of   a   satellite  is 
proportional  to   the  cube  of  its   mean  distance  from  the 
center  of  the  earth. 

THE    ELLIPTICAL    ORBIT 


The   size  and   shape   of   the   elliptical   orbit  is   deter- 
mined  completely   by   the  parameters   a   and   e.    (See   Fig. 
2-1).      These  quantities   can  be   expressed  in  terms   of 
other   orbital   dimensions   when  a  and   e  are  not   given 
directly.      For   example,    the   ellipse  is   determined  if    only 
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Figure  2-1.    The   Elliptical   Orbit 
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the  minimum  and  maximum  altitudes   of  the  satellite  are 
known.      The  perigee   distance  p  is   given  by  the  minimum 

altitude  h    .      plus   the  radius   of   the  earth   R.      Similarly. 

min  F  J  * 

the  apogee  distance  q  is   given  by  h  +   R.      The  semi- 

max 

major  diameter  a  is   then   (p+q)/2  and  the   eccentricity  e, 

the   semi-minor   diameter  b   and   semi-latus   rectum   L    are 

given  respectively  by 

e   =  ^lE  (2-1) 

a 


b  =  aVl-e2  (2-2) 

L  =  -  (2-3) 

a 

from   the  geometric   properties   of   the  ellipse.      The  angular 

displacement    &    of  the  satellite  from  the  perigee  is 

measured  in  the  direction  of   satellite  motion  and  is 

commonly  called  the  true  anomaly. 

It  remains   to   determine  the  position  of   the  satellite 

in  its   orbit  as  a  function  of  time.      In  polar  coordinates, 

this   position  is   given  by 

r    -  a(l-e2)  (2-4) 

1+e  cos'Q 

Unfortunately   Q     is   not  a   simple  function  of  time, 

and  the   satellite  position  is  more  easily   expressed  in 


terms   of   the   eccentric   anomaly    E,    whose  geometric 
relationship   to   the   satellite  position  is   illustrated  in  Fig. 
2-2. 

From    Newton's   laws,    the  radial  acceleration 

(9.80665    m/sec    )    of  a  body   in   circular   orbit   of   radius   R 

2 
is    equal   to   v     /R.       Solving  for   v     yields   a  velocity  of 

7,90i+»455   meters/second. 

Then  the  period  in  seconds  is   given  by 

p^   =  £7/r_  =   510^.6536   seconds.  (2-5) 

o  v 

vo 

By  Kepler's  third  law,    the  period  of  any  other  near   earth 
satellite  is  then 

P   =   P      r^l  seconds,  (2-6) 

°\1 

and  the  mean  angular   motion  ^1  is  given  by 

p.  =  -p—     radians/second.  (2-7) 

Now   define  the  mean  anomaly  as  follows, 

M  =  u   (t-T)  (2-8) 

where  t  is  time  in  seconds  and  T   is  the  time  the  satellite 
passes  through  perigee.      From   this,    Kepler's   equation 

M  =    E  -e  sin  E  (2-9) 

can  be  solved  for   E.    (See   Chapter  III).      From   Fig.    2-2, 
it  is  apparent  that 

r  =  a    (1  -e  cos   E)  (2-10) 
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2  2 

=    SD^    +   DC 


'!  S'd]2    +    (OC-OD)2 

—  a  sin  El      +    (ae-a   cos   E) 
I 
=  a2    (1-e2)    sin2   E  +  a2    ( e-cos   E)2 

=  a2    (1-e  cos   E)2 


Figure  2-2.    Position  vector  in  terms   of   E 
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from    which,    together   with   equation    (2-i+) 

\ 


cos 


e  = 


cos    E   -e 
i  -e   cos    E 


sin 


e  =  «vi  -e2 


sin   E 


(2-11) 


Equations    (2-11)    make  it   possible  to    write  the   cartesian 
coordinates   of   S   in  Fig.    2-2   as 


-I    = 


cos  G    =  a    (cos   E  -e) 
e?n  6     =   a  VI   -e        sin   E, 


r   sin 


(2-12) 


GEOCENTRIC  EQUATORIAL  COORDINATES 

The   reference   system   to   which   the   orbital   parameters 
are  referred  is  known  as   the  geocentric   equatorial 
coordinate  system.      It  is  a  conventional  right-handed 
cartesian  system   with  the  origin  at  the   center  of  the 
earth.       (See  Fig.    2-3).      The  x-axis   is   directed   toward 
a  fixed  point  in  space  where  the  sun's   apparent  path 
through  the  heavens   crosses   the  earth's   equatorial  plane 
in  a  northerly  direction.      This   point  is  variously  known 
as   the  first   point   of  Aries   or   the  position  of   the  vernal 
equinox.      It   is   not   truly  fixed  in  space  as    stated  in  the 
basic   assumptions,    but   its   movement   due  to   nutation  of 
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tellite 


HT  direction  of  vernal   equinox 

XL  longitude  of  ascending  node 

CJ  argument   of  perigee 

i  inclination  of   orbit 

Q  true  anomaly 


Figure  2-3.    Projection  of   Orbit   on  Unit 
Geocentric   Sphere 
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the   earth's   spin  axis   is   only  a   small  fraction  of   a   degree 
per  year  and   can  be  ignored  in  this   problem.      The  y-axis 
is  90°    from  the  x-axis,    measured  eastward  along  the 
equator,    and  the  z-axis  is  in  the  direction  of   the  north 
pole.      Corresponding  unit   coordinate  vectors  are  designated 

The  orientation  of  the  orbit  plane  in  the  reference 
system   is  given  by  the  parameters  fl  ,  cJ   ,    and  i.      The 
longitude  of  the  ascending  node  fl  is   the  angle   measured 
eastward  along  the  equator  from  the  x-axis   to  the  point 
where  the  satellite  crosses   the  equator  plane  in  a  northerly 
direction.      The  argument   of  the  perigee  O  is  the  angle 
measured  in  the  orbit  plane  in  the  direction  of   motion 
from   the  ascending  node  to  the  perigee  position.      The 
inclination  i  is  the  angular  separation  of  the  equator  and 
orbit  planes.      This  angle  is  less   than  90°   if  the  satellite 
motion  is   eastward    (or  direct).      If  the  motion  is   west- 
ward  (retrograde)    the  angle  is   the  supplement  of  the 
angular  separation  of  the  planes. 

The  change  from  x1,    y'   coordinates   in  the  orbital 
plane  to  x,y,z   coordinates   can  be  accomplished  by  the 
linear  transformation 

12 


x  =    Px   x»   +    Qx  y» 


y  =    Py    x'   +   Qy  y» 
z   -   Pz   x'   +    Qz   y' 


(2-13) 


where  the  P's  and   Q's  are  the  direction  cosines   of  the  x1 
and  y1  axes   relative  to  the  x-axis    I2J  .      These  six 
quantities  are  found  to  be 


P  =  cos  CJ     cos  JT    -   sin  CJ     sin  -0-    cos  i 

Qx  =  -sin  U)     cos-A-  -  cos  CJ     sin  -TL    cos  i 

P  =  cos  oJ     sinJl     +  sin  D      cos  S\    cos  i 

Q  =  -sin  cJ     sin.0-     +  cos  oJ     cos  -TL  cos  i 

P  =  sin  u)     sin  i 
z 

Q  =  cos  u)    sin  i. 


\ 


>•      (2-U) 


The  eccentric  anomaly   may  be  introduced  directly 
into   equations    (2-13).      If   we  define 


Ax  =  aPx>         B: 

Av   =  aPv,         B. 


=  aVTU 


Az   =  aPz,        Bz   =  aVl   -e2      Qz 


=  aVTTe2 
=  aV7Te^ 


Q- 


Q, 


then 


x  =  A  (cos    E  -e)    +   B      sin   E 

y  =  A  ( cos   E  -e)    +  B      sin  E 
j  y 

z   =  A,  (cos    E   -e)    +   B      sin   E 

Z  z 


(2-15) 


(2-16) 
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The  first  time  derivatives   of  these  position 
coordinates  yield  the  components  of  the  satellite  velocity 
vector  in  terms   of   E.      They  are  found  to  be 

~    =    4Sp    (-Av   sin   E   +   Bv   cos    E)      ) 
dt  r  x  x 

&2L    =    JS&   (_A„  sin  E  +  BTr  cos    E)       }'  (2-17) 

at  v  j  y 

<||-    -    M.    (-Az   sin   E   +   Bz   cos    E)      j 

CONSERVATION   OF    ENERGY,    MOMENTUM  AND 
AREAL   VELOCITY 

Kepler's   second  law   and  the  laws   of   conservation  of 

energy  and  angular  momentum  reveal  relationships   which 

will  be  found  useful.      The  famous  Vis  Viva  integral,    the 

energy   equation  for  an  elliptical  orbit,    states   that 

v2   =   GM     - i-|     meters2/second2  (2-18) 

where  G  is  a  gravitational   constant  and  M  is  the  mass  of 
the  earth.      In  the  MKS  system   of   measurement,    the 
product  GM  is   equal  to  3.98  x   10        meters-y seconds     and 
if  nautical   miles  are  used  as  units   of  distance,    GM  is 
equal  to     6.267  x   ICk      n.m.    /seconds    .      It  is   noticed 
that  for  any  orbit  with  a  given  major  axis,    v  is  a  function 
only  of  r,    or   equivalantly,    of  altitude. 
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Angular   momentum   is   conserved  in   the   two   body- 
problem,    and   its   vector   representation  is    given  by   the 
constant   vector 

H    =   r   x   —  .  (2-19) 

dt 

Multiplication  by   r«  yields   the   equation   of   the   orbital  plane r 

r    •    H   -  O.  (2-20) 

From    Kepler's    second  law,    the  areal   velocity   of   a   body  in 
an   elliptical   orbit   is    constant    I2I  .      It   can  be   shown  that 

^    =1h.  (2-21) 

dt  2 

Then  by  taking   the  time   integral   of   both   sides,    the  area 
swept   out  in  time   t   is   given  by 

A=   i- Ht   +    constant.  (2-22) 

2 

In  particular,    if   the  area   swept   out   between   two   points   r1 
and   r?   is   known,    than  the  time  in  the   orbit  between  these 

two   points   is 

2 
r21        ff  a21- 
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CHAPTER   III 
SATELLITE    ORBIT    COMPUTATION 

SELECTION   OF    ORBITAL    ELEMENTS 

A  satellite  orbit   may  be  determined  from  a   wide 
variety   of   given  parameters.      In  the  present   model, 
those  initial  conditions  and  parameters   are  set   which  the 
analyst   may  consider  to  be  of  particular  importance.      Of 
first   interest   are  the   maximum   and   minimum   altitudes, 
from   which  the  size  and   shape  of  the  orbit  are  determined. 
Specifying  the  inclination  immediately  determines  the 
northern  and  southern  limits   of  latitude  over   which  the 
satellite  may  pass  and  the  direction  of   motion.      Setting 
the  argument  of   the  perigee  fixes  the  perigee  point  in 
the  desired  hemisphere,    north  or  south.      Finally,    the 
time  at  perigee  and  longitude  of   the  ascending  node  are 
chosen  or   computed  as  indicated  below  and  the  orbit  is 
computed  from   the  equations   of   Chapter  I. 

POSITION  VECTORS  IN  TERMS  OF  TERRESTIAL 
LATITUDE  AND  LONGITUDE 

The  analyst   may  wish  to  know  the  components   of  a 

position  vector   corresponding  to  a  position  on  earth 
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designated  by   latitude  and  longitude  at   a   given  time.      For 
example,    he   may   wish  to   specify  an  initial  geographic 
position  for   the  perigee  of  an   orbit   or   determine  the 
position  vector   of    the  interceptor   launch   site.      Conversely, 
he   may   wish  to   know   the   geographic   location  of   a   satellite 
position  when  its   position  is   given  as   a   vector. 

The  Greenwich  Hour  Angle   of   the  first   point   of  Aries 
(GHAT    )    is   listed  in  the   Nautical  Almanac   for   the  begin- 
ning  of   each   day   of   the  year    [lj  •      GHAf      is   defined  as 
the   central  angle   measured   westward  along  the   equator 
from   the   Greenwich   meridian  to   the  vernal    equinox. 
Designate  the  longitude   of   a   position  by  \      with   eastern 
longitudes   being  positive,    and   the  latitude  by    0     with 
northern  latitudes   being  positive.      Then  the   Sidereal   Hour 
Angle    (SHA)    between   the  vernal   equinox  and  the   position's 
meridian  is    equal   to    GHAT    +  X  .      However    SHA  increases 
at   the  rate  of   the   earth's   rotation,    which  in  the   chosen 
coordinate   system   is    equal   to    one  revolution  per   sidereal 
day  or    .0000729   radians/second.      If  the  value  of   GHAHQ 
is   taken  as   that   given  for   O      January   1   of   the  year,    and 
t  is   the  number   of   seconds   after   the  beginning   of   the  year, 

SHA  =   GHAX   +A   +   0.000729t    radians  (3-D 
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and 

x  =   r   cos   SHA   cos    o 

y  =   r   sin    SHA   cos   S        >«  (3-2) 

z    =   r    sin  6 

If  the  position  vector  r(x,y,z)  is  given  the  process 
may  be  reversed  although  care  must  be  taken  in  assigning 
signs   to   the   computed  angles. 

r     =    yx2+y2+z2  (3-3) 

&    =      sin"1  If)  (3-4) 

SHA  =   sin^^^J  (3-5) 

•where   -   -L—   ^    £  ^    -2—    with   the   same   sign  as    z   and 
-Tf    ^      SHA   ^  if     with   its    magnitude  and   sign   determined 
by  seeing  which  quadrant   contains   the  x,y  coordinates   of 
r.        A  may  then  be  found  directly  from   equation   (3-1). 

If,  as  suggested  above,  the  analyst  wishes  to  specify 
an  initial  perigee  point  and  longitude  of  ascending  node  from 
geographical  coordinates,  the  corresponding  position  vectors 
r".  and  i^_  are  determined  from  equations  (3-1)  and  (3-2). 
XI    is   of    course   simply   equal   to    SHA  and 

<*)    =   cos"1     ILdiiJ  (3_6) 

r      r 

with   the   same   sign  as    z.      It    must   be  noted  that   these 
two  points  fix  the  orbital  plane  so  that  the  inclination   may 
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no   longer   be   chosen  arbitrarily.    However,    it   is    easily 
computed,    for   the   vector 

*I  =  5vx  F<o  (3-7) 

is  inclined  to   the   z-axis   by  the  angle  i   if    the   satellite 

If 

motion  is   direct   and        — -  +   i   if   the   motion  is   retrograde. 


Then 


1    (or  j+   i)    =   cos  r>       i 

SOLUTION    OF    KEPLER'S    EQUATION 

Given  the   orbital  parameters  and   the  satellite's 
position  in  the  orbital   plane,    the   computation  of   coordinates 
in  the  geocentric    equatorial   coordinate   system   is   a   straight- 
forward  process.      However,    the  position  of   the   satellite 
in  its   orbital   plane  as   a  function  of   time  is   found   from 
the   solution  of    Kepler's    equation  which   is   transcendental 
in  E.      This    equation  is   solved  by  an  iterative   method   of 
successive  approximations. 

Write   equation    (2-9)    as 

E   =   M   +   e   sin   E^  (3-8) 

o  ■ 

and   setting   E      equal   to   M,    solve  for    E.       Next   setting   E 
equal   to   the   computed  value   of    E,    repeat   the   computation. 
The  procedure  is   repeated  until   the   desired  accuracy  is 
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obtained.      It   can  be  shown  that   the  iterated  values   of 
E   converge  rapidly   to   a   single  value. 
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CHAPTER  IV 
INTERCEPTOR  TRAJECTORY  COMPUTATION 

The  first  problem   to  be  considered  may  be  stated  as 
follows.      If  an  interceptor  launch  position  and  a  predicted 
intercept  point  are  given,    how   can  one  determine  trajectory- 
parameters   that   will   yield  a   successful  interception.      In 
particular,    out   of  all  possible  interceptor  trajectories, 
which   one  requires   the  least   energy  to   be   supplied  to   the 
interceptor  vehicle. 

After  the  interceptor  has   exhausted  its  propulsive 
energy,    it  is   subject   only  to  the  same   central  force  of 
gravitation  as  is  the  satellite.      It  will  follow  an  'elliptical 
orbit   with  a  major   semi-diameter  determined  by  its  altitude 
and  velocity  in  accordance  with   equation   (2-18).      Conversely, 
if  the  semi-diameter   of  the  interceptor  orbit  is  known, 
the   energy  of   the   orbit   is   fixed.      Thus   the   orbit   with 
the  minimum   semi-diameter   which  passes   through   the  launch 
site  and  the  intercept  position  is   the   minimum    energy 
interceptor   trajectory. 

This   reasoning  assumes   that  the  interceptor  attains 
its   maximum   velocity  immediately  after  its   launch  and 
thereafter  follows   a  purely  ballistic  path   with  air  resistance 
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neglected.      Actually  there  is  an  initial  phase  of  its  flight 
during   which  it   acquires   the  velocity  and  attitude 
necessary  for   the  proper   ballistic   path   to   the  intercept 
position.       This   phase   depends   upon  the  thrust   character- 
istics  of  the  propulsion  system,    the  payload  weight  and 
other  factors   of  interest   to  the  designer  and  engineer. 
This   paper   considers   only   the  ballistic   requirements   which 
the  propulsion   system    must   meet,    and   the  interceptor   is 
taken  to   be  injected   into   its   orbit   with   the  necessary 
velocity  vector  at   the  time  of  launch. 

Consider   the  plane   determined  by  the  position  vector 
of  the  launch  site  and  the  position  vector  of  the  intercept 
point   as   shown  in  figure  4-1. 

Let   C  be  the   midpoint   of   line   BI    and   let   the 
distance   CF=    1/2    GI .       (Notation:    BI    represents   the 
magnitude  of   the  directed   line   segment   BI.)       Note  that 
by   construction  the   distance   OG+GI+IF   =   OB+BF,       Now 
define   e^  as   the  ratio   Cl/CF   and  p^   as   the   quantity 
CF(eh2   -1).      Also   let   bh2   =   CF(ph).       Construct   the 
hyperbola   HH'    from   the   equation 

rh  =  _Ph t*rx) 

1   -   eh   cos  Q 

22 
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center   of  the  earth 


launch   site 


I  intercept  point 

R         radius   of   the  earth 

h  satellite  altitude 


Figure  4-1  •    Locus   of   Secondary  Foci  of 
Interceptor   Trajectories 
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•with  the  origin  at  I .      B   is  the  focus   of  a   mirror  hyper- 
bola and  from  the  geometric  properties   of  the  hyperbola, 
the  distances  from   I   and  B   to  any  point  on  the  curve 
( F'  for   example)    differ  by  a  constant,    in  this   case  equal 
to   GIo      Thus  the  distance  OB+BF'  is   equal  to  the 
distance  OI+IF'.      But  recall  that  the  sum   of  the  distances 
from  any  point  on  an  ellipse  to  the  two   foci  is   a  constant, 
equal  to   the   major  diameter  2a.      Therefore  any  allipse 
with   one  focus  at  the  center  of  the  earth  and  the  other 
lying  on  the  hyperbola  HH1   will  pass  through  the  launch 
site  and  the  point  I .      Clearly  the   ellipse  with  the  smallest 
value  of  a  is   the   one  where   F=F!,    thus   determining  the 
minimum   energy  orbit. 

The  orbital   elements   of  the  minimum   energy  inter- 
ceptor trajectory   may  now  be  found  from  the  given 
position  vectors   OB  and  OI  •      BI=OI   -  OB.      The  distance 
GI   is   the  altitude  h  of  the  satellite  at  intercept  and  is 
equal  to   OI   -  R. 

OF    =    OB    +    BF  (4-3) 

a   =    OB    +    BF  (4-4) 
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OF  ,,     -» 

e   =  -5 •  (4-5) 

2a 


The  vector   product   OK  =   OB   x  OI    is   normal   to   the  plane 
of   the   orbit,    and  its  component    will   be  positive  if   the 

motion  of   the  interceptor   is   direct.      Then 


-1       OK    .    k    \  ,,     ,, 

i   =   cos  \.  (4-6) 

1         OK       I 

The  vector   product   ON  =   k  x   OK  points   in  the   direction 
of  the  ascending  node,    so 

It  =  cos"1    (ON    :    Li  (4-7) 

I        ON        I 

Perigee   distance  is   equal   to   a  -   OF,    so   the   perigee  position 
vector 

OF   =    (OF   -a)    OF  (4-8) 

and  the  final   spatial   element 


0)   =   cos"1      ,ONw°Pv   •  U-9> 

The  period  of  the   orbit   and  the   mean   motion    p.    of 
the  interceptor  are  found  from   equations    (2-6)    and   (2-7). 
The   sine   of   the   true  anomaly   of   the  interceptor   at   the 
launch   site  is   given  by 
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.     q  OP    x   OB  ,,     ... 

sin  U     =     , — -'  (4-10) 

(OP)(OB) 

and  from    equation    (2-11),    sin   E   and   E   are   determined. 
Equation    (2-9)    gives   the  value  of   M   and  finally 

t    -T    =21  (4-H) 

M 

from    which   the   epoch   of   the   orbit   is   determined.      A 

similar   computation   will   give  the  time  from   perigee  to 

intercept   from    which   the  time   of   flight   of   the  interceptor 

is   known. 

The  above   derivation   determines   the   orbital   parameters 
for   the   minimum    energy   orbit  but   the   same  procedure   may 
be   followed   whenever   the   location   of   the   secondary  focus 
is   given.      The   orbital   plane  is   determined  by  the   location 
of   the  launch   site  and   the  intercept   point,    and  any   secondary 
focus   position  in   that   plane   can  be   computed  by   equation 
(4-1)  •      A   transformation  of    coordinates    will  yield  the 
position  vector   of   F'   from    which   the   desired   parameters 
may  be   computed  by  the   methods   above. 

The  velocity  vector   can   now   be   computed  by   equation 
(2-17)    for   any  point   in  the  orbit.      Of    course,    the  points 
of   greatest   interest  are   the  point   of   interception  and  the 
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launch    site.      At   the  intercept   point,    the   difference 
between  the   satellite  and   interceptor   velocity   vectors   is 
the   closing  velocity.      The   velocity  vector   at   the   launch 
site  is   the   sum   of   the   velocities   imparted  by   the  propulsion 
system   and   the  velocity  vector   of   the  launch   site  due  to 
the   earth's   rotation.      At   lower   latitudes,    this   latter 
velocity  approaches    1000    mph   and  an  interceptor   with   a 
given  propulsion   system    will  have  a   substantially  longer 
maximum   range  when  launched  to   the   east   than   when 
launched   toward  the   west.       The  velocity  vector   of   the 
launch   site  is    seen  from    equation    (3-2)    to   be 

—  =   -uR    sin   SHA   cos  6 
dt 

77=  uR   cos   SHA   cos   &  (4-12) 

dt 

dz 

dt   -   ° 

This   quantity   must  be  subtracted  from   the   orbital   velocity 
vector  at   that   point   if   the  velocity   to   be   supplied  by  the 
propulsion   system    alone  is   to   be  determined. 

Figure  4-2    shows   the   minimum    energy   orbit   for   the 
conditions   displayed   in   Figure  4-1   together   with   another 
interception  orbit  having   a   secondary  focus    F'    on  the 
hyperbola    HH'. 
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Figure  J+-2.    The  Minimum    Energy   Trajectory 
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CHAPTER  V 
COMPUTER  PROGRAMMING  CONSIDERATIONS 

The   satellite  interception  problem,    with  its    complex 
and   constantly   changing   relationships   of   the  many   para- 
meters  and  variables   involved  lends   itself   very   well   to 
simulation  by  a   modern  high   speed  computer.      Huge  amounts 
of  information  can  be   computed,    stored  and  recalled  with 
great   speed  and  accuracy.      Such  a  computer   model  has 
been  constructed  for   execution  on  a   Control   Data   Corp- 
oration computer   model  1604*      (See  Appendix  A). 

Substantial   economies   can  be  effected  by  writing 
computer  programs   in  machine  language  in  terms   of 
execution  time  and  memory  requirements.      This   comes  at 
the  expense  of  using  a  format  which  is   not   easily  adapted 
for  use  by  other   computers  and  is  not  readily   modified  by 
other  programmers.      An  alternative,    which  is   the  approach 
used  in  this   model,  is  to   write  the  computer  program   in  a 
commonly  used  compiler  language,    a  variation  of   FORTRAN 
known  as  FORTRAN-63    [J  .      The   model  follows   the 
mathematical  approach  outlined  in  this  paper,    and  by 
noting  the  liberal  number   of   explanatory  comments   in  the 
program,    an  analyst   should  be  able  to  adapt  the  program 
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for  his   own  purposes   on  another  computer   employing  a 
FORTRAN    compiler   language.      The  remainder   of   this 
chapter  is   devoted  to  inspection  of   computer  and   compiler 
language  considerations   which  affect  the  construction  of 
the  model. 

Inasmuch  as   the  satellite  orbit  is  fixed  in  the 
chosen   coordinate  system,    much  computational  time  can 
be  saved  by  computing  position  and  velocity  vectors  for 
one  orbital  period  and  storing  the  information  for  recall 
as   needed.      If  sufficient   storage  space  is  available  in  the 
computer  memory,    all  quantities   of  interest  such  as 
position  vector   components,    velocity  vector  components, 
true  and  eqcentric  anomalies   with  their  trigonometric 
functions,    distance  from  the  center  of  the   earth,    etc., 
may  be  stored  for  later  use,    and  this  is   recommended 
where  possible.      However,    the  period  of   a   circular   orbit 
at   1000    miles  altitude  is  nearly  7500   seconds,    and  storage 
of  all  quantities   of  interest  in  sufficiently  small  increments 
of   time  imposes  a   considerable  requirement  for   memory 
space. 

In  this   particular   model,    only  the  position  vector 
components   can  be  stored  for   each  second  of  time  in 

30 


normal  format    without   exhausting   the   memory   storage 
capacity.      By   going  outside   the   FORTRAN-63   language 
to   a   compatible  assembly  program    CODAP-1,    it   is 
possible  to   store  the  eccentric  anomaly  for   each   second 
by  storing  two   quantities  in  a  single  storage  address    I5J  » 
The   CODAP-1   subroutine   BYTE5   which  accomplishes   this 
is  not  included  in  the  FORTRAN-63   listing. 

A  small  bank  of  temporary  variable  names    (CI   through 
C9)    in  a  COMMON  statement   saves   some  storage  by  using 
the  same  names  wherever  an  intermediate,    non-permanent 
quantity  is   computed.      The  same  statement   contains 
conversion  quantities  that  are  needed  in  several  sub- 
routines.     Additional  information  may  be  stored  on  magnetic 
tape  or   equivalant  devices  but  this  is   profitable  only  if  the 
time  required  for  access   to  the  information  is  less   than 
the  time  required  to   compute  it. 

Liberal  use  is   made  of   subroutines  and  functions  in 
the  interest  of  flexibility,    for  the  model  is   constructed 
so  that   more  relationships  than  are  explicitly  specified  may 
be  considered  with  only  minor   modifications.      For   example, 
the  model   can  easily  consider  a  terminal  guidance  aspect 
of  the  problem   or  serve  as  the  basis  of  a  probabalistic 
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approach  to   interception.       Communication  between  the 
subroutines   is   facilitated  by  the  use  of   COMMON    state- 
ments. 

Clarity  of  the  program   for  other  analysts  has  been 
enhanced  by  consistent  use  of   suffixes  and  prefixes  in 
variable  names  and  the  use  of   mnemonic  variable  names 
throughout. 
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appendix  a 
fOrtran-63  computer.. simulation  program 

The  program    listed  here   has   been   extensively   tested 
for  a   wide  variety   of   orbital    elements   and  the  generated 
orbits   show    close  agreement   with   the   orbits   of  actual 
satellites.      Storage  is   provided  for   satellite   positions  in  one 
second   increments   for   orbital  altitudes   to    1000    miles. 
Sample   output    statements   have  been  included   to    suggest 
ways   in   which   the  user   might   call   for   computed  values. 

The   computational   procedure   follows   the   methods 
outlined  in  the   thesis   with   extensive  use   of   subroutines    and 
functions   for   those   calculations   that   are  used  repeatedly 
or   in  other   sections   of   the   program .      This    enhances   the 
flexibility   of   the  program   for   the  analyst   who    may   want 
the   calculated   quantities   for   supplementary   considerations 
in  the  interception  problem. 
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